#include <stdio.h>
#include <math.h>
#include "wiz_cephes.h"

static const double rel_error = 1E-12;

double MACHEP = 1.11022302462515654042E-16;  // 2**-53
double MAXLOG = 7.09782712893383996732224E2; // log(MAXNUM)
double MAXNUM = 1.7976931348623158E308;      // 2**1024*(1-MACHEP)
double PI     = 3.14159265358979323846;      // pi, duh!

static double big    = 4.503599627370496e15;
static double biginv = 2.22044604925031308085e-16;

int sgngam = 0;

double wiz_cephes_igamc(double a, double x)
{
    double ans, ax, c, yc, r, t, y, z;
    double pk, pkm1, pkm2, qk, qkm1, qkm2;

    if ((x <= 0) || (a <= 0))
        return (1.0);

    if ((x < 1.0) || (x < a))
        return (1.e0 - wiz_cephes_igam(a, x));

    ax = a * log(x) - x - wiz_cephes_lgam(a);

    if (ax < -MAXLOG)
    {
        return 0.0;
    }
    ax = exp(ax);

    /* continued fraction */
    y    = 1.0 - a;
    z    = x + y + 1.0;
    c    = 0.0;
    pkm2 = 1.0;
    qkm2 = x;
    pkm1 = x + 1.0;
    qkm1 = z * x;
    ans  = pkm1 / qkm1;

    do {
        c  += 1.0;
        y  += 1.0;
        z  += 2.0;
        yc  = y * c;
        pk  = pkm1 * z - pkm2 * yc;
        qk  = qkm1 * z - qkm2 * yc;
        if (qk != 0)
        {
            r   = pk / qk;
            t   = fabs((ans - r) / r);
            ans = r;
        }
        else
            t = 1.0;
        pkm2 = pkm1;
        pkm1 = pk;
        qkm2 = qkm1;
        qkm1 = qk;
        if (fabs(pk) > big)
        {
            pkm2 *= biginv;
            pkm1 *= biginv;
            qkm2 *= biginv;
            qkm1 *= biginv;
        }
    } while (t > MACHEP);

    return ans * ax;
}

double wiz_cephes_igam(double a, double x)
{
    double ans, ax, c, r;

    if ((x <= 0) || (a <= 0))
        return 0.0;

    if ((x > 1.0) && (x > a))
        return 1.e0 - wiz_cephes_igamc(a, x);

    /* Compute  x**a * exp(-x) / gamma(a)  */
    ax = a * log(x) - x - wiz_cephes_lgam(a);
    if (ax < -MAXLOG)
    {
        return 0.0;
    }
    ax = exp(ax);

    /* power series */
    r   = a;
    c   = 1.0;
    ans = 1.0;

    do {
        r   += 1.0;
        c   *= x / r;
        ans += c;
    } while (c / ans > MACHEP);

    return ans * ax / a;
}


/* A[]: Stirling's formula expansion of log gamma
 * B[], C[]: log gamma function between 2 and 3
 */
static unsigned short A[] = {
    0x6661, 0x2733, 0x9850, 0x3f4a,
    0xe943, 0xb580, 0x7fbd, 0xbf43,
    0x5ebb, 0x20dc, 0x019f, 0x3f4a,
    0xa5a1, 0x16b0, 0xc16c, 0xbf66,
    0x554b, 0x5555, 0x5555, 0x3fb5};
static unsigned short B[] = {
    0x6761, 0x8ff3, 0x8901, 0xc095,
    0xb93e, 0x355b, 0xf234, 0xc0e2,
    0x89e5, 0xf890, 0x3d73, 0xc114,
    0xdb51, 0xf994, 0xbc82, 0xc131,
    0xf20b, 0x0219, 0x4589, 0xc13a,
    0x055e, 0x5418, 0x0c67, 0xc12a};
static unsigned short C[] = {
    /*0x0000,0x0000,0x0000,0x3ff0,*/
    0x12b2, 0x1cf3, 0xfd0d, 0xc075,
    0xd757, 0x7b89, 0xaa0d, 0xc0d0,
    0x4c9b, 0xb974, 0xeb84, 0xc10a,
    0x0043, 0x7195, 0x6286, 0xc131,
    0xf34c, 0x892f, 0x5255, 0xc143,
    0xe14a, 0x6a11, 0xce4b, 0xc13e};

#define MAXLGM 2.556348e305


/* Logarithm of gamma function */
double
wiz_cephes_lgam(double x)
{
    double p, q, u, w, z;
    int    i;

    sgngam = 1;

    if (x < -34.0)
    {
        q = -x;
        w = wiz_cephes_lgam(q); /* note this modifies sgngam! */
        p = floor(q);
        if (p == q)
        {
        lgsing:
            goto loverf;
        }
        i = (int)p;
        if ((i & 1) == 0)
            sgngam = -1;
        else
            sgngam = 1;
        z = q - p;
        if (z > 0.5)
        {
            p += 1.0;
            z  = p - q;
        }
        z = q * sin(PI * z);
        if (z == 0.0)
            goto lgsing;
        /*      z = log(PI) - log( z ) - w;*/
        z = log(PI) - log(z) - w;
        return z;
    }

    if (x < 13.0)
    {
        z = 1.0;
        p = 0.0;
        u = x;
        while (u >= 3.0)
        {
            p -= 1.0;
            u  = x + p;
            z *= u;
        }
        while (u < 2.0)
        {
            if (u == 0.0)
                goto lgsing;
            z /= u;
            p += 1.0;
            u  = x + p;
        }
        if (z < 0.0)
        {
            sgngam = -1;
            z      = -z;
        }
        else
            sgngam = 1;
        if (u == 2.0)
            return (log(z));
        p -= 2.0;
        x  = x + p;
        p  = x * wiz_cephes_polevl(x, (double *)B, 5) / wiz_cephes_p1evl(x, (double *)C, 6);

        return log(z) + p;
    }

    if (x > MAXLGM)
    {
    loverf:
        return sgngam * MAXNUM;
    }

    q = (x - 0.5) * log(x) - x + log(sqrt(2 * PI));
    if (x > 1.0e8)
        return q;

    p = 1.0 / (x * x);
    if (x >= 1000.0)
        q += ((7.9365079365079365079365e-4 * p
               - 2.7777777777777777777778e-3)
                  * p
              + 0.0833333333333333333333)
             / x;
    else
        q += wiz_cephes_polevl(p, (double *)A, 4) / x;

    return q;
}

double
wiz_cephes_polevl(double x, double *coef, int N)
{
    double  ans;
    int     i;
    double *p;

    p   = coef;
    ans = *p++;
    i   = N;

    do
        ans = ans * x + *p++;
    while (--i);

    return ans;
}

double
wiz_cephes_p1evl(double x, double *coef, int N)
{
    double  ans;
    double *p;
    int     i;

    p   = coef;
    ans = x + *p++;
    i   = N - 1;

    do
        ans = ans * x + *p++;
    while (--i);

    return ans;
}

double
wiz_cephes_erf(double x)
{
    static const double two_sqrtpi = 1.128379167095512574;
    double              sum = x, term = x, xsqr = x * x;
    int                 j = 1;

    if (fabs(x) > 2.2)
        return 1.0 - wiz_cephes_erfc(x);

    do {
        term *= xsqr / j;
        sum  -= term / (2 * j + 1);
        j++;
        term *= xsqr / j;
        sum  += term / (2 * j + 1);
        j++;
    } while (fabs(term) / sum > rel_error);

    return two_sqrtpi * sum;
}

double
wiz_cephes_erfc(double x)
{
    static const double one_sqrtpi = 0.564189583547756287;
    double              a = 1, b = x, c = x, d = x * x + 0.5;
    double              q1, q2 = b / d, n = 1.0, t;

    if (fabs(x) < 2.2)
        return 1.0 - wiz_cephes_erf(x);
    if (x < 0)
        return 2.0 - wiz_cephes_erfc(-x);

    do {
        t   = a * n + b * x;
        a   = b;
        b   = t;
        t   = c * n + d * x;
        c   = d;
        d   = t;
        n  += 0.5;
        q1  = q2;
        q2  = b / d;
    } while (fabs(q1 - q2) / q2 > rel_error);

    return one_sqrtpi * exp(-x * x) * q2;
}


double
wiz_cephes_normal(double x)
{
    double arg, result, sqrt2 = 1.414213562373095048801688724209698078569672;

    if (x > 0)
    {
        arg    = x / sqrt2;
        result = 0.5 * (1 + erf(arg));
    }
    else
    {
        arg    = -x / sqrt2;
        result = 0.5 * (1 - erf(arg));
    }

    return (result);
}
